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Abstract. Initially introduced in the framework of quantum control, the so-called monotonic 
algorithms have demonstrated excellent numerical results when dealing with various bilinear 
optimal control problems. This paper presents a unified formulation that can be applied to 
more nonlinear settings. In this framework, we show that the well-posedness of the general 
algorithm is related to a nonlinear evolution equation. We prove the existence of the solution 
to this equation and give important properties of the optimal control functional. Finally we 
show how the algorithm works for selected models from the literature. 



1. Introduction 

This paper aims at presenting a general unified formulation of several algorithms that were 
proposed in different areas of nonlinear (bilinear) control. Given a cost functional J{v) depending 
on the control v, these algorithms are iterative procedures that construct a sequence of solution 
candidates v'^ with the important "monotonic" behavior, i.e. J{v''^^) < J{v'^) ; the algorithms 
have been named after this property as " monotonic" . An convenient property of these procedures 
is that the monotonicity does not requires any additional computational effort, but results from 
the construction of the procedure itself. 

These procedures have first been used in the field of quantum control, where quantum particles' 
dynamics is controlled by a laser field as described by the Schrodinger equation (cf. Section UTT] 
for more details about the modeling of this problem). In this framework, the function that 
to each control associates the final state of the system is highly nonlinear. This induces poor 
performance of standard, gradient-based algorithms. The "monotonic schemes" introduced in [3J 
I49[ 155] appeared in this context and were found to perform excellently in this very nonlinear 
setting. These schemes were used in bi-linear situations where the control multiplies the state. 
These were soon followed variants [54l[56] that generalized the cost functional to include situations 
more complicated that a distance to a given target. 

At first the relationship between the procedures introduced in the cited works was not obvious 
but in [29] it was showed that there are all particular cases of a two-parameter class of algorithms. 

Though the monotonic schemes are based on algebraic calculations, the specific setting induced 
by the Schrodinger equation enables in [43] to relate the monotonic schemes to trajectory tracking 
algorithm [30, STj. At the numerical level, efficient discretizations of the procedure have been 
proposed in (281 ITT] and a time parallelized version was introduced in [37] . Continuing interest in 
the monotonic schemes lead to the introduction in [T3] of versions that ensure that the resulting 
field will fit in a given frequency window. 

In previous works the objective was encoded through a criterion on the final state but adapta- 
tions were proposed in [36] |35] to deal with the case where the optimal control functional depends 
on the whole dynamics of the control process or when the dynamics involve integro-differential 
equations. Additional situations involving dissipation operators were proposed in [32j with a 
non-Markovian version in 1331 . 
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Further different examples consider the case where the system is described by a density matrix 
operator instead of a wave function: details on the computation and convergence proofs limited 
to this case were given in [331 1371 |3H] and in [TH] this is applied to create a quantum computer 
gate; further applications can be found in [37l |9]. 

At some point similar procedures were also proposed in other control applications f[8l 119). 

The convergence of the algorithm has been obtained in the case of quantum control (see 
Section 14. 1|) using Lojasiewicz-Simon inequality (see [71 [171 EH SS] and the references therein) 
and also in discrete and continuous settings in [31 [31] . The structure of the proofs shows that 
when J is analytic and its gradient is Fredholm, convergence is guaranteed as soon as J contains 
a penalization term of the L^-norm of the control v. Note also that another proof has been 
presented in the framework of semi- group theory jl5l using compactness arguments. All these 
results are available for a bi-linear setting when the control multiplies the state and are specific 
to the Schrodinger equations. 

Up to this point all works presented above considered bi-linear situations (in all cases the opti- 
mal control functional is a non-linear functional) ; only recently different cases were documented 
in the literature: in [42l [20] the procedures were tailored to tackle specific non-bilinear models 
in which the control field appears up to power 3. A situation when a unique control appears at 
arbitrary powers of polynomial was proposed in [34] . A model where the system is a nonlinear 
Bose-Einstein condensate was given in |46| . 

In all situations where monotonic algorithms were introduced the well-posedness of the algo- 
rithms were proved by ad-hoc techniques and the same for convergence, although the algebraic 
computations share similar points. The purpose of this paper is to identify the similarities present 
in all these situations, and present a general setting to which the "monotonic" algorithms belong ; 
we also propose corresponding formulas and procedures to solve such type of problems. This al- 
lows to tackle general non-linear situations the cannot be solved with techniques presented in the 
literature. 

The paper is structured as follows: Section [2] defines the general framework where our pro- 
cedure applies. The algorithm itself is presented in Section |3l At this point we show that the 
well-posedness of the algorithm is related to a nonlinear evolution equation, and prove the exis- 
tence of the solution to this equation. We also give important properties of the optimal control 
functional. Some examples of concrete realizations follow in Section [4] together with numerical 
results illustrating the application of the algorithm. 



2. Problem formulation 

Let E, T-L and V be Hilbert spaces with V densely included in %. We denote by -e and (•, •)v 
the scalar product associated with E and V. 

For any two vector spaces A and B we denote by C{A, B) the space of linear continuous of 
operator between A and B. 

Given a real or complex valued function we denote by V x'f its gradient with respect to 
the variable x. We also denote by and D^^x the first and the second derivative of vectorial 
functions in the Frechet sense. 

Remark 1. Recall that, given Hi and H2 two Banach spaces and U C Hi an open subset of Hi, 
a function f : U H2 is said to be Frechet differentiable at x ^ U if there exists a continuous 
linear operator £ ^{Hi, H2) such that 

y \\f{x + h)-fix)-Axih)\\H, 

hm —■ = 0. 

h^O \\h\\H, 
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// this is the case the operator Ax is called the Frechet differential ( or Frechet derivative ) of f at 
X and is denoted Ax ~ F>xf. 

Let us also recall that for a open set f2 C and any Hilbert space Hi, L°°(f2; Hi) is the space 
of functions f from with values in the Hilbert space Hi such that for almost all t € ft the norm 
\\f(t)\\Hi is bounded by the same constant (the lowest of which will be the L°°{Q; Hi) norm of f). 

In the same way one can define L^{fl;Hi) 

(1) L^{n;Hi) = {f -.n^ Hi such that [ \\f{t)\\%^dt < 00} . 

When the derivatives of functions of f are considered the Sobolev spaces W^'°° have to be 
introduced; we refer to [TJ I53j for further details. 

Within an optimal control formulation, the evolution of a system X(t) is encoded in the 
following optimization problem: 

(2) minJ(w), 

V 

where 

(3) J{v):^ ( F{t,v{t),X{t))dt + G{X{T)). 

Jo 

The functions F:Mxi?xV— >]R and G : V — > M are assumed to be differentiable and integral 
assumed to exist. The system is described by a state function X{t) E V being governed by the 
evolution equation 

(4) dtX + A{t,vit))X = B(t,v{t)) 

(5) X{Q) - Xo. 

where v : [0,r] _E is the control. The unbounded operator A{t,v) : W x E x H ^ H is such 
that for almost all t € [0, T] the domain of A{t, v)^^^ includes V; furthermore we take B{t, v) such 
that for almost all t e [0, T] and all w € £; we have B{t, v) G £('H, n£(V, V* ). We postpone to 
Section[3] (cf . Lemma [3^ Theorem[T]) the precise formulation of additional regularity assumptions 
to be imposed on A, B, F, G. 

Remark 2. Finally, note that E is not necessarily a real number, neither finite dimensional, cf. 
Section \4.2\ This means that the control can be a set of several time- dependent function but also 
a distributed control depending on time and also on a spacial variable. 

Let us stress that although the equation is linear in X (for v fixed) the mapping w i— > X is not 
linear ; the term A{t,v{t)) multiplies the state X and as such the mapping is highly nonlinear 
(of non-commuting exponential type). 

Remark 3. Most of the previous works considered a bilinear operator A{t,v) i.e., A(t,v)X = vX; 
the only exceptions (cf. discussion in the Introduction) were of the polynomial type (of order at 
most 3 in \42\ I20j and polynomial with E = M.^ in [34j ). The techniques present in the above 
papers cannot be used for general operators A{t,v). On the contrary the results in this paper 
include all the situations considered in the bibliography but also apply to all nonlinearities in v 
compatible with the hypothesis of Lemma \3.4\ and Thm. [7] below. 

Moreover, the following concavity with respect to X will be assumed throughout the paper: 

(6) VX, X' e V, G{X') - G{X) < (VxG(X), X' - X)v, 



(7) Vt e M, Vi) e E, yX, X' e V, F{t, v, X') - F{t, V, X) < {VxF{t, v, X),X' - X) v 
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Remark 4. Contrary to the more technical hypothesis that will he assumed latter, the properties 
([7]) and the linearity of (j4]) are crucial to the existence of the monotonic algorithms. 

3. Monotonic algorithms 

We now present the structure of our optimization procedure together with the general algo- 
rithm. 

3.1. Tools for monotonic algorithms. The monotonic algorithms exploit a specific factoriza- 
tion which is the consequence of the results in this section. To ease the notations we will make 
explicit the dependence of X on u, i.e. we will write X„ instead of X in Eqs. (|4H5]). 
We define the adjoint state (see [131 [21]) by: 

(8) dtY,- A*{t,v{t))Y,+VxF{t,v{t),X,{t)) = 

(9) Y,{T)=VxG{X,{T)). 
A first estimate about the variations in J can be obtained: 
Lemma 3.1. For any v',v : [0,r] — > E denote 

T(t,X,{t),vit),v'{t),Y,{t),X,,{t)'^ ^ -{Y,{t), (^A{t,v\t)) - A{t,v{t))^X,,{t))^ 

(10) {Y,{t),B{t,v'it)) - B{t,v{t)))v + F{t,v'it),X„,it)) - F{t,vit),X^,{t)). 
Then 

(11) J{v')-J{v)< T{t,X,it),vit),v'{t),Y,it),X,,it)y 
Proof. Using successively (HI),©, (HJ and finally ([9]), we find that: 

Jiv')-J{v) = [ F{t,v{t),X,'it)) - F{t,vit),X4t)) 



]dt. 



+ F{t, v'{t), X,, [t]) - F{t, v{t), X,, {t))dt 
+G{X^,{T)) -G{X^(T)) 

< [ {VxF{t,vit),X,it)),X,,it)-X,it))^ 
Jo 

+ F{t, v'{t),X.,, (t)) - F{t, v{t), X,: {t))dt 



+ {Y,{T),X,,{T)~X,{T)) 



V 



< 



{jYvit) - A{t,v{t)yY,{t)+S/xF{t,v{t),X,{t)),X,,{t) - X,{t))v 



V 



- {Y,{t),(^A{t,v'{t))-A{t,v{t)))x,,{t)) 
+ {Y,{t),B{t,v'{t))-B{t,v{t))U 

+ F{t,v'{t),X,,{t)) ~ F{t,v{t),X,,{t))dt. 

Due to ([8|), the first term of the left-hand side of this last inequality cancels and the result 
follows. □ 
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Remark 5. The focus of the result is not on obtaining an estimation of the increment J{v') — J{v) 
via the adjoint (which is well documented in optimal control theory, cf. |13[ I25j ); we rather 
emphasis that the evaluation of the integrand T(...) at time t requires information on the control 
v{s) for all s G [0,r] (in order to compute X^{T) then Y^{t)) hut on the second control v'{s) 
only for s G [0,t] (because this is enough to compute Xyi{t)). This estimate can be useful in 
deciding, at time t if the current value of the control v' (t) will result in an increase or decrease 
of J{v'). This localization property is a consequence of the concavity of F and G (in X) and 
bi-linearity induced by A. The purpose of the paper is to construct and theoretically support a 
general numerical algorithm that exploits this remark. 

Remark 6. We can intuitively note that T(...) has the factorized form: 

(12) T(t, Xy{t), v{t), v'{t),Y,{t), X,, (t)) = A{v, v'm -e {v'{t) - v{t)) , 

with ■ E the E scalar product. Thus v' can always be chosen so as to make it negative ( in the worse 
case set it null by the choice v' — v). We will come back with a formal definition of A{v,v'){t) 
and a proof of the previous relation in Section \3.3\ 

A more general formulation can be obtained if we suppose that the backward propagation 
of the adjoint state is performed with intermediate field v (cf. also [29] ), i.e. according to the 
equation : 

^Yj;-A*{t,v{t))Yv + VxF{t,v{t),Xy{t)) = 
Y^iT)^VxG{X,iT)). 

Note that because of its final condition, Yy actually also depends on v. Nevertheless, for sake of 
simplicity, we keep the previous notation. We then obtain the following lemma. 

Lemma 3.2. For any v',v,v : [0,T] — > E, 

J{v')-J{v) < ~{Yr,{t),(A{t,v'{t)) - A{t,v{t)))Xy,{t)U 

+ {Yy{t),B{t,v'{t)) - B{t,v{t)))^ 
+ F{t,v'{t),Xy,{t)) ~ F{t,v{t),Xy,{t))dt 
+ -{Yyit),(^A{t,vit)) ~ A{t,vit))^Xyit))^ 
+ {Yy^t) (t), B{t, vit)) - B{t, v{t))U 

+ F{t, v{t), Xy{t)) - F{t, v{t),Xy{t))dt. 

In this lemma, the variation in the cost functional J is expressed as the sum of two terms, and 
can be considered as factorized with respect to v' — v and v — v. 

Remark 7. Lemmas \3.1\ and \3.2\ are generalization of previous results that were proved in the bi- 
linear case. To the best of our knowledge, only specific corollaries requiring additional assumptions 
have been used in the literature up to now. 
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3.2. The algorithms. The factorization obtained above enables to design various ways to ensure 
that J(w') < J(v), i.e. that guaranty the monotonicity resulting from the update v' ^ v. This 
allows to present a general structure for our class of optimization algorithms. We focus on the 
one that results from Lemma |3. II 



Algorithm 1. ( Monotonic algorithm ) 

Given an initial control , the sequence {v^)k£n is computed iteratively by: 

(1) Compute the solution X^k of with v — v'' . 

(2) Compute the solution Y^k of with v = v'' , starting from 

Y,k{T) ■.= VxG{X,k{T)). 

(3) Define v^'^^ together with X^k+i such that for all t < T the following monotonicity con- 
dition be satisfied : 

(13) Aiv''+\v'')it) -E (v^+^t) ~ v\t)^ <0. 

Lemma lOI then guarantees that J(u'^+^) < J{v^). Several strategies can be used to ensure (fT3|): 
we will present one below. Its importance stems from the fact that no further optimization is 
necessary once this condition is fulfilled. In order to guarantee many authors (see [^1^155] ) 
consider an update formula of the form: 

(14) v''+\t)-v''{t)^-^A{v''+\v''){t), 

where is a positive number, that can also depend on k and t. In what follows, we focus on the 
existence of solution of p^ . and on practical methods to compute it. If v''^^ satisfies p^ . the 
variations in J satisfy: 

j(v'^+i) - j(f'=) < ~e r {v^+\t)-v^{t)fdt. 

Jo 

Note that p4p reads as an update formula combining on the one hand a gradient method: 

„fe+l(^)_„fe(^) = _lA(^;^^;'=)(^), 

and on the other hand the so-called Proximal Algorithm (introduced by [6 ), which prescribes: 

v''+Ht)-vHt) = -]-A{v^+\v''+^){t). 


Remark 8. When F = and A is independent of v, i.e. linear control with final objective, (|14p 
coincides with a gradient method. 

3.3. Well-posedness of the algorithm. In this section, we focus on the procedure obtained 
when using Algorithm[l]with the update formula (|14|) . Since this procedure involves the resolution 
of an implicit equation, see Eq. (fM)) , we prove the existence of a solution and present a convergent 
procedure to compute it. As a by-product, we obtain a proof of the monotonicity of the algorithm. 

Lemma 3.3. Suppose that for any t S [0,T]; 

-A-.R-xYxVxE-^R defined by A{t,X,Y,v) = {Y, A{t,v)X)yr is Frechet differentiahle 
everywhere with respect to v for any X,Y,v. 

-S:MxVx_E— >R with B{t,Y,v) — {Y, B{t,v))\- is Frechet differentiate everywhere with 
respect to v for any Y, v. 

- F is Frechet differentiable everywhere with respect to v £ E for any X,Y,v. 
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Then there exists A(-, ■■,t,X,Y) e C°{E'^,E) such that, for all v,v' e E 

A{v',v;t,X,Y)-E {v' -v^ = - (y, (^Ait,v') - Ait,v)'^X + Bit,v') - B{t,v)) 



V 



(15) +F{t,v',X)-F{t,v,X). 

Moreover, if A,B,F are of class in v then IS.{-,-]t,X,Y) can he defined through the explicit 
formula: 



A{v',v-t,X,Y) = ^ -V^(^{Y,A{t,w)X ~ Bit^w))^) 
(16) +\I^F{t,v + \{v' -v),X)d\. 



w—i)-\-X{v' — v) 



(17) /^^{v',v) := ' {v' -v)eE. 



Proof. We denote by || • || the norm associated with E. Since A,B^F are Frechet difTer- 
entiable with respect to v the fuU expression in Ec^. (jl5p is of the form S(u') — S(w) with 
S(u) = —A{t, X, Y, v) + B{t, Y, v) ~ F{t, V, X) differentiable in v; we introduce 

\\v' — v\\ 

Since S is differentiable, we obtain the continuity of Ah(w', v) for aU points v' = v and A.s{v, v) 
V„S(ti) (the continuity is obvious everywhere else) hence the conclusion. 
Finally, Eq. (|16p is an application of the identity 

E{v')-E{v)= [ \/,E{v + X{v' -v))dX-E{v' -v). 



□ 

Lemma 3.4. Suppose that 

- A,B,F are of (Frechet) class with respect to v with DwA, DwB uniformly bounded as 
soon as X ,Y are in a bounded set; 

- 'S/vF is of class in X ; 

- DyyF{t, ■, X) is bounded by a positive, continuous, increasing, bounded from below function 

X <^ H\\x\\). 

Given e > 0, {t,v,X,Y) eRx E xV xY and a bounded neighborhood W of (t, v, X, Y), there 
exists 9* > depending only on e, W , \\v\\, \\X\\ and \\Y\\ such that, for any 9 > 9* 

(1) A(t)', v; t, X, Y) = —9{v' — v) has an unique solution v' = Vg{t, v, X, Y) G E. 

(2) Vg{t,v,X,Y) implies 

(18) - V, ( (y, A{t, v)X)^ ) (v) + V„ ( (y, Bit, «)) V ) (v) + VyF{t, V, X) = 0. 

(3) \\Ve{t,v,X,Y)-v\\ < "^llll'^ll+IIJII+fe(ll^ll) {Afo(t) + Mi||t;||} with Mo{t) and Mi indepen- 
dent ofv,X,Y. If the dependence of A,B,F on t is smooth then Mo{t) is bounded on 
[0,T]. 

(4) ^dit,v, X,Y) is continuous on W . 

(5) Let X belong to a bounded set; then X i— >■ Ve{t,v,X,Y) is Lipschitz with the Lipschitz 
constant smaller than e. 

Proof. 

(1) Denote h — v' — v and Qt,v,x,Y{h) = ~^('"+^^'''*'^'^) _ When the dependence is clear we 
will write simply G{h) instead of Gt,v,x,YW- We look thus for a solution to the following 
fixed point problem: Q{h) — h. For 9 large enough, the mapping C/ is a (strict) contraction 
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and we obtain the conclusion by a Picard iteration. The uniqueness is a consequence of 
the contractivity of Q. 

(2) Ev' = v then h = thus g{h) =0 which gives (HH) after using ([TO)) . 

(3) For 9 large enough, the mapping Q is not only a contraction but has its Lipschitz constant 
less than, say, 1/2. Because of the contractivity of t?, we have — ||tj(0)|| < 1?(0)|| = 
\\g{h) - g{Q)\ < which amounts to < 2||t?(0)||. Next, we note that 

ii^.nMi <- II A(t;, V, t, X, Y) - A(0, 0, t, X, Y)\\ + || A(0, 0, t, X, Y)\\ 

\\y[Oj\\ < 7 



< Ah\\v\\ + Ahit) 

and the estimates follows. 

(4) Formula shows that A depends continuously on t,v,v',X,Y. Consider converging 
sequences i„ — t, v„ — >■ v, X„ — )■ X, F„ — >■ F and define /i„ := Vg{tn,Vn, Xn,Yn) and 
h := Vg{t,v,X, Y). 

Given W and ?7 > 0, consider large value of 9 such that: 

- for any {t' ,v' , X' ,Y') £ W, Qti,v',X',Y' is a contraction with Lipschitz constant less 
than 1/2. 

- for any {f ,v' ,X' ,Y'), {t" ,v" , X" ,Y") e W, 

II A(«" + h, t", X", r") - A(z;' + ft, t', X', Y')\\ < 77. 

This last property implies \\Gt„,v„,x„,Y„{h) ~ Gt,v,x,Y{h)\\ < ^ for n large enough. On 
the other hand 

\\hn-h\\ = \\Gt„,v„.X^,Y^{K) - Gt,v,X.Y{h)\\ 

< \\Gt„,v„,X„,Y„{hn) — Gt„,v„,X„,Y,Ah)\\ 
+ \\Gt„,v„.X„,Y„{h) - Gt,v.X,Y{h)\\ 

< l\\hn~h\\ + ^. 

We have thus obtained that for n large enough : ^||/i,i — /i|| < ^ and the continuity 
follows. 

(5) Subtracting the two equalities 

A{Vi,v;t,Xi,Y) = -9{Vi-v), A{V2,v;t, X2,Y) = -9(y2 - v) 

and using that A(y, v; t, X, Y) is in X and v gives to first order 

Av{...){Vi - V2) + Ax{...){Xi - X2) = -9iVi - V2). 

For 9 large enough the operator Ay(...) -\- 9 ■ Id\s invertible and the conclusion follows. 

□ 

Remark 9. Note that 9* is proportional to (||X|| v|l>^|| v + ll^^llv + fc(||X||v)). 

We are thus able to give an example of a setting where the existence of v'^^^ (t) satisfying ((Ti 
is guaranteed. 



Theorem 1. Suppose thatA,B,F satisfy hypothesis of Lemma \S.4\ Also suppose that the opera- 
tors A, B are such that Eqs. and have solutions for any v £ L°°{0, T; E) with v 1— >■ X , 
V 1-^ Y locally Lipschitz. 
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(1) For any v £ L°°(0,T; E), there exists 0* > such that for any 9 > 9*, the (nonlinear) 
equation 

(19) dtX,' (t) + A{t, v')X,, (t) = Bit, v') 

(20) v'{t)^Vg{t,v{t),X,,{t),Y4t)) 

(21) X,,{0)^Xo 

has a solution. Here Yy is the adjoint state defined by flH^i and corresponding to control 

V. 

(2) There exists a sequence {9k)keN such that the algorithm (cf Section[ 

a/ initialization v'^ G L°°{0,T; E), 
b/ v'^+^it) = Vs, (i, w'=(t), (t), Yy. [t)) 

is monotonic and satisfies 



V 



L2([o,r])- 



(3) With the notations above, if for all t £ [0,T] v'^'^^{t) = ^'^(t) (i.e. algorithm stops) then 
v'^ is a critical point of J: V„J(i''=) = 0. 

Proof. Most of the proof is already contained in the previous lemmas. The part that still has 
to be proved is the existence of a solution to (IT9l) - ((2T|) . 

Given v G L°° {0,T; E), consider the following iterative procedure : 

vo = V, vi+i{t) = Vg{t,v{t),Xy^{t),Yy{t)). 

We take a spherical neighborhood By (R) of v of radius R and suppose 

Vfc <l, VkE By{R). 

Since the correspondence v i— > Xy is continuous, it follows that the set of solutions Sy,R := 
{Xyy-jW G By{R)} of ^ IS boundcd. In particular for w = vi hy the item [3] of Lemma [3.41 the 
quantity \\Ve{t,v{t), Xy^{t),Yy{t)) — v\\ will be bounded by R for 9 large enough (depending on 
R, independent of /), i.e. vi+i S By{R). Thus vi S By{R) for all I > 1. 

Since Sy^n is bounded, recall that by itemlSlof Lemma l3.4l the mapping X i-> Vo{t, v{t), X, Yy(t)) 
has on Sy^R a Lipschitz constant as small as desired. Since the mapping w i— >■ X^ is Lipschitz, 
for 9 large enough, w G By{R) i— ^ Ve{t,v(t), Xyj,Yy{t)) is a contraction. By a Picard argument 
the sequence vi is converging. The limit will be a solution of ([T9l - [20l) . □ 



4. Examples 

We now present two examples that fit into the setting of Theorem [TJ The space does not allow 
to treat all other variants (cf. references in Introduction) so we leave them as an exercise to the 
reader. 

Within the framework of control theory, nonlinear formulations prove useful nowadays in do- 
mains as diverse as the laser control of quantum phenomena (see [24l [38l [39l |40l [5TJ [52]) or the 
modeling of a equilibrium (or again social beliefs, product prices, etc) within a game with infinite 
numbers of agents (see [U] [551 [53]). Yet other applications arise from modern formulations of 
the Monge-Kantorovich mass transfer problem (see [5] [4] [8]). 

In the following, we present some examples coming from these fields of application and present 
the corresponding monotonic algorithm resulting from Theorem [1] 



4.1. (I): Quantum control. 
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4.1.1. Setting. The evolution of a quantum system is described by the Schrodinger equation 

dtX + iH{t)X ^ 
X(0,z)=Xo(z), 

where i = V--T, H{t) is the Hamiltonian of the system and z G M? the set of internal degrees 
of freedom. We assume that the Hamiltonian is an auto-adjoint operator over L'^ {M."^ ; C) , i.e. 
Hit)* = Hitj^. Note that this results in the following norm conservation property 

(22) 11-'^ (^7 •)IU2(R^:C) = I|-'^o||l2(rt.c), yt > 0, 

so that the state (or wave-) function X(t,-), evolves on the (complex) unit sphere 5* := |x e -L^(R'^; C) : H-'f Ili2(][{ 

The Hamiltonian is composed of two parts: a free evolution Hamiltonian Ho and a part that 
describes the coupling of the system with an external laser source of intensity v{t) £ 'R, t > 0; 
a first order approximation leads to adding a time-independent dipole moment operator /i(x) 
resulting in the formula H(t) = Hq — v{t)^i and the dynamics: 

dtX + i{Ho-v{t)^)X = Q 
X(0) = Xo. 

The purpose of control may be formulated as to drive the system from its initial state Xq 
to a final state Xtarget compatible with predefined requirements. Here, the control is the laser 
intensity v{t). Because the control is multiplying the state, this formulation is called "bilinear" 
control. The dependence v H> X{T) is of course not linear. 

The optimal control approach can be implemented by introducing a cost functional. The 
following functionals are often considered: 

(23) J{V) := \\X{T) - Xtarget\\h(^tL...C) + / a{ty{t)dt, 

Jo 

(24) J{v) := -(X(T),0(X(T)))i.(K.^c) + C a{ty{t)dt, 

Jo 

where O is a positive linear operator defined on Ji, characterizing an observable quantity and 
a{t) > is a parameter that penalizes large (in the sense) controls. The goal is thus to 
minimize these functionals with respect to v. According to ()22|) the cost functional J is equal to 

(25) Jiv):=2-2Re{X{T),Xtarget)LHR-r:C)+ f a{ty{t)dt, 

Jo 

so that the functionals J and J satisfy assumptions ([6]) and 0. 

4.1.2. Mathematical formulation. We have 

• A{t^v) = Ho + v(t)fi with (possibly) unbounded u-independent operator Hq (but which 
generates a semi group) and bounded operator /i. The dependence of A on v is smooth 
(linear) and therefore all hypotheses on A are satisfied. 

• E ^ R, n = L2(]R'';C), V = dom{Hl'^) (over C), or their realifications V. = x L^, 
V = dom{Hl'^) X dom{Hl'^) (over M) see [16]; 

• B{t,v) = 0. 

• F{t,v,X) — a{t)v{t)'^ with a{t) £ i°°(M); here the second derivative D^yF is obviously 
bounded. Since it is independent of X it will be trivially concave. 



For any operator M, we denote by M* its adjoint. 
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• G is either (see, e.g., [211 [29]) 2 - 2Re{Xtarget, X{T))yr or -{X{T),OX{T))yr where O is a 
positive semi-definite operator; both are concave in X . 

• Here 



(26) 



A{v', V] t, X, Y) = ~Re{Y, ifiX)^ + a{t)(v' + v) 



and the equation in v' lS.(v' ,v]t, X,Y) = —6{v' — v) has, for 6 large enough, an unique solution 
v' = Vg{t, V, X, Y) (s-<^it))v+Iie{Y,^^.xw 

• at the k + 1-th iteration. Theorem [T] guarantees the existence of the solution X'^^^ of the 
following nonlinear evolution equation: 

(27) idtX^+\t) = (^Ho 

Then 

(28) 



k+i(,. frr , {0-a{t))v'+Re{Y,.,tfiX'^+')^ , 
fc+i _{e- a{t))v'' + Re{Y^, , i^lX''+^}v 



a{t) 



4.1.3. Numerical test. In order to test the performance of the algorithm we have chosen a case 
already treated in the literature. The system under consideration is the O — H bond that vibrates 
in a Morse type potential V{x) — £)o(exp(— /3(x — x')) — 1)^ — Dq. The dipole moment operator 
of this system is modeled by /i(x) = ^o.xe~^ We refer to [55l for more details concerning this 
system. The objective is to localize the wavefunction in a time T = 131000 at a given location xq ; 
this is expressed through the requirement that the functional J is maximized, where the observable 
O is defined by 0{x) = .^g-ioi^-^^oy _ y^i^ consider a constant penalization parameter a = 1 

and optimization parameter = 10^^. The numerical values we use are give in the next table. 



Do 


13 


x' 


X* 


Xo 


70 




0.1994 


1.189 


1.821 


0.6 


2.5 


25 


3.088 



Results are presented in Fig. [T] 



" iteraltions 



20 25 



Figure 1 . For the example in Section 14.11 we plot the evolution of the cost 
functional values J{v^) as function of the iteration number k. Monotonic decrease 
is observed as expected by the theoretical arguments. 



4.2. (II) : Mean field games. 
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4.2.1. Setting. Although the Nash equihbrium m game theory has been mitiaUy formulated for a 
finite number of players, modern results (see [211 [22l [23] ) indicate that it is possible to extend it to 
a infinite number of players and obtain the equations that describe this equilibrium; applications 
have already been proposed in economic theory and other are expected in the behavior of multi- 
agents ensembles and decision theory. 

The equations describe evolution of the density X{t, z) of players at time t and position z € 
Q = [0, 1] in terms of a control v{t, z) and a fixed parameter i> > 0: 

dtX ~ vAX + div{v{t, z)X) = 0, 
X{0) = Xo. 

The control v is chosen to minimize the cost criterion For reasons related to economic 
modeling interesting examples include situations where F, G are concave in X, e.g., as in jl9| 

(29) G = 0, F{t, z,X)^ f pm - Pz)X{t, z) + 1^ ' + ^^(i, z)dz, 

Jq Ci+C2X{t,z) 2 

with positive constants /3, co,ci,C2 and p(i) a positive function. Another example is given in [8]: 

(30) G{X{T))= / V{z)X{T,z)dz, F{t,z,X)= / X{t, z)v^{t, z)dz, 



where V encodes a potential. The interpretation of this terminal cost is that the crowd aims at 
reaching zones of low potential V at the terminal time T. 

The relevance of the monotonic algorithms to this setting has been established in several 
works [8lfT9]. 

4.2.2. Mathematical formulation. We have 

, E = W^^°^{0, 1), H = L^{0, 1), V = H\0, 1) see [l9] and [10] (Chap XVIII §4.4) 

• A{t,v) = —i^A ■ +div{v-). The dependence of A on is smooth (linear) and therefore all 
hypotheses on A are satisfied {Dy^A = 0, ...). 

• B{t,v) = 0. 

• with definitions in 1^ F{t,v,X) = jQp{t){l - l3z)X{t,z) + + ^^^^X{t, z)dz; 
F is concave in X; the second differential DyyF has all required properties. 

• G — (algorithm will apply in general when G is concave with respect to X). 

• Here 

(31) A{v',v;t,X,Y)=VY+^^ 

and the equation in v' A{v' ,v;t, X,Y) = —9{v' — v) has for all > an unique solution v' = 
Vg{t,v,X,Y) := — . 

• at the k + 1-th iteration. Theorem [T] guarantees the existence of the solution X^^^ of the 
following nonlinear evolution equation: 

(32) dtX'^+\t) - vAX''+^ + dz^;( tlZ?l!^-_Z^x'=+i) = 0. 

+ 1/2 

Then 

(33) V - g ^ , A„fe+i-A 

4.2.3. Numerical test. The algorithm is test on the time interval [0, 1] with p{t) = 1 and the 
numerical values /3 = 0.8, cq — C2 — 1 , ci = 0.1. Results are presented in Fig. [2] 
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Figure 2. For the example in Section [4.21 we plot the evolution of the cost 
functional values J{v^) as function of the iteration number k. Monotonic decrease 
is observed as expected by the theoretical arguments. 



4.3. Additional application. As a third example we consider a nonlinear vectorial case from [T^l 
[50] which differs from that of Section O in that v{t) = (^^^ e E = R'^ and A{t, v) = 

i[Ho + {viit)^ + V2it)^)ni + vi{t)^V2it)fi2]- Here, denoting = ~Re{Y,iiJ,iX)v + a{t), ^2 = 
—Re{Y,ifi2X)-\/ we obtain 

(34) A«,,x.y)=6(-t:;)+6(<-,;5y") 



v' = Ve{t,v,X,Y) 



and the equation in v': A(w', v; t, X, Y) = —6{v' — v) has, for 9 large enough, an unique solution 

/ (6>-Ci)t>2-6f? \ 

, ) 

?-Ci+C2 '"^''V^l ^^"^ • We leave as an exercise to the reader the 

\ ^1^1^ (<'-ei)i-2-S2"^ / 
\ 0+51+6 eqrj^^ / 

writing of the equation for X^^"^ and the formula for 

This model corresponds to the problem of controling the orientation 7 of a molecule, considered 
as rigid rotator. 

4.3.1. Numerical test. To test our approach we have used the parameters of the molecule CO [T^ 
[50] . namely i/o — BJ^, where B is the rotational constant and J is the angular momentum. We 
consider the basis given by the spherical harmonics ; the corresponding matrix is diagonal with 
diagonal coefhcients given by (i?o)fc,fc — k{k + 1). The controlled is performed over an interval of 
lenght T = 20Tper — 20-g. We consider constant penalization factor a = 10~^ and optimization 
parameter 9 — 10"^. 

The other parameters correspond to the polarizability and the hyperpolarizability components of 
the molecule. In this way we have ni = — ^A, and fi2 = — f /3, with A — ^(A|| cos^ 7 + Aj^ sin^ 7), 
(3 — — 3/3_l) cos'^ 7 + 3/3_L cos 7). The matrix cos 7 is tridiagonal, with: 

(cos7)fe,fc = 0, 

fc + 1 



(cos7)fc,fc+i = (cos 7) 



k+l,k 



^(2/c + l)(2fc + 3) 
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We use the numerical values given in [12l [50] : 



B 




^11 


/3,| 


P± 


1.93 


11.73 


15.65 


28.35 


6.64 



The results are presented in Fig. [3] As in the previous examples, a rapid convergence is 
obtained, since about 100 iterations are necessary to reach the numerical convergence. 




20 40 60 .m 100. 120 140 160 180 200 

iterations 



Figure 3. For the example in Section 14.31 we plot the evolution of the cost 
functional values J{v^) as function of the iteration number k. Monotonic decrease 
is observed as expected by the theoretical arguments. 



5. Conclusion 

Motivated by a set of control algorithms that were initially introduced in the specific context 
of quantum control we have presented an abstract formulation that includes them all. It is seen 
that the algorithm involves at each step a highly nonlinear evolution equation. We identified the 
theoretical assumptions that ensure that the evolution equation is well posed and has a solution. 
The proof being constructive it serves as basis for numerical approximations of the solution. 
We also proved several properties concerning the algorithms and more specifically concerning 
its convergence. Examples are provided to indicate how the procedure proposed solves previous 
cases from the literature and also new situations that were not previously considered. Numerical 
simulations indicate that the procedures have indeed the expected behavior. 
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